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We study the switching-process of the magnetization in a ferromagnetic-normal-metal multilayer 
system by a spin polarized electrical current via the spin transfer torque. We use a spin drift- 
diffusion equation (SDDE) and the Landau-Lifshitz-Gilbert equation (LLGE) to capture the coupled 
dynamics of the spin density and the magnetization dynamic of the heterostructure. Deriving a fully 
analytic solution of the stationary SDDE we obtain an accurate, robust, and fast self-consistent 
model for the spin-distribution and spin transfer torque inside general ferromagnetic/normal metal 
heterostructures. Using optimal control theory we explore the switching and back-switching process 
of the analyzer magnetization in a seven-layer system. Starting from a Gaussian, we identify a 
unified current pulse profile which accomplishes both processes within a specified switching time. 

PACS numbers: 75.76.+J, 72.25.Pn, 75.70.Cn, 85.75.-d 



I. INTRODUCTION 



Spin transfer torque in nanoscaled ferromagnetic/normal-metal (FN) heterostructures has potential application 
for data storage and manipulation [JH5]. Apart from the experimental studies many theoretical investigations have 
been made since the pioneering work by Slonczewski and Berger [3HS]- The problem to describe the physics in 
FN heterostructures arises from the need to consider the dynamics of the conduction electrons as spin carriers and 
the dynamics of the localized magnetic moments in parallel and in different regions of the heterostructure. The 
electron dynamics is faster by several orders of magnitude than that of the latter [S]. Moreover, the spin dynamics 
in normal metal regions differs significantly from that in ferromagnetic regions: the former is characterized by fast 
diffusion and slow spin relaxation, while in the latter the opposite is the case. This time hierarchies make it difficult 
to provide a fully numerical solution. A Boltzmann-transport theory for magnetic multilayer systems including 
the spin was developed by Valet-Fert [101 111] . On the next level of approximation a drift-diffusion equation was 
applied for mobile spins [12]. The dynamics of the localized magnetic moments is governed by the Landau-Lifshitz 
Gilbert equation (LLGE), extended by additional spin transfer terms. A similar investigation has been performed for 
semiconductor/ferromagnetic multilayers assuming ballistic transport, but using non-equilibrium Green's functions 

In this paper we utilize this time hierarchy and base our model on an exact stationary solution to the spin drift- 
diffusion equation which we were able to obtain for constant electric current and arbitrary but piece- wise constant layer 
parameters. We solve self-consistently the LLGE and the spin drift -diffusion equation (SDDE) for the conduction 
electrons in an external magnetic field to explore switching scenarios as a function of current pulse profiles. 

The paper is organized as follows. In Sec.[n]we present our model for FN multilayer system. Sec. |III| and Sec. |IV| are 
devoted to the mathematical description of the magnetization dynamics (LLGE) and the dynamics of the conduction 
electrons (SDDE) respectively. The exact solution of our SDDE is given, with details deferred to the Appendix. In 
Sec. [V] we present numerical results for a symmetric seven-layer system. Optimal current pulse profiles to switch the 
magnetization in a given time from parallel to antiparallel state (and in opposite direction) is shown. Our results are 
compared with our fully numerical simulations to confirm the validity of our approach. Our results regarding critical 
switching currents versus switching time agree well with earlier work by others [HI 115) . 
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II. MODEL 



Our model of the heterostructure assumes three different physical building blocks: (i) the normal-metal leads and 
spacer layers, (ii) ferromagnetic polarizers , and (iii) ferromagnetic analyzers. The leads and the spacer layers are 
chosen to be nonmagnetic (N) metals with equal material properties. A lead is assumed to be infinitely thick and 
serving as a spin bath with vanishing spin polarization. We describe a wide ferromagnetic hard polarizer layer (Pi, 
P2 in Fig. [I]) as static and homogeneous. A thin ferromagnetic (soft) analyzer layer (region A in Fig. [IJ is treated 
as a ferromagnetic mono-domain described by a single time-dependent variable, a unit-vector m(t) pointing in the 
direction of the magnetization [THl Q~7]. The conduction spin-electrons are treated as classical magnetic moments 
moving in an external magnetic field created by localized magnetic moments in the fcrromagnet. The spin density 
S(x, t) is the dynamical variable to describe the spin distribution [18] . It is defined for an isolated ferromagnet with 
magnetization direction m as 

S = nV~m. (1) 

Here ri = 11^ + is the free electron number density, where n^^ is the particle density with spin up/down respectively 
and V = corresponds to the spin density polarization, extracted from the experiment |19| . In this work we use 

for the spin density the dimensionless quantity s = S/nS. For simplicity we do not consider spin-resolved quantities 
but use mean values instead (diffusion constant, electric conductivity, spin diffusion length etc.). 



III. MAGNETIZATION DYNAMICS 



A. Landau-Lifshitz-Gilbert equation 



The temporal evolution of the magnetization M is governed by the LLGE [TBI [27Jj • Using the saturation magneti- 

' lm\ - _J_ (9M\ 
at ) s t ~ m s \ dt ) st 



zation M s , we define the quantities M = M s m, h = 7H, where 7 is the gyromagnetic ratio and ( — — 1 



to obtain an equation of motion for the dimensionless magnetization: 



dm 1 , a , , s f <9m\ 

, m x h -m x (m x h) + — — . (2) 



dt 1 + a 2 1 + a 2 ' \ dt 

Here h = h an + h ex is the effective field containing the anisotropy field and external fields measured in units of a 
frequency and a the Gilbert damping constant. With a unit vector n we set for the anisotropy field 



uj an n(m ■ n) , (3) 
= £m x (AI S x m) . (4) 



where u> an is the corresponding frequency, (^f) denotes the spin-transfer term [41 [51 117). 



9m 

~dt 

Here AI S = (I s ) i„ — (I s ) out stands for the spin current absorbed inside the domain, whereas £ is a constant [2"T1 122j . 
Without external torque the equilibrium magnetization is either parallel (P) or antiparallel (AP) to n. 



B. Dipole field 

In this paper we consider the control of the magnetization by spin currents only. So the only contribution to h ex 
from the outside are the dipole fields originating from the polarizers. In order to obtain a simple result and a crude 
estimate of the order of magnitude of the dipole fields we consider a polarizer (here written for Pi in Fig. [I]) as a 
cylinder with radius R and thickness X\ which is homogeneous magnetized and compute the field at the position 
x m = (x% +Xs)/2. Evaluation of the general integral for a dipole density |23j we obtain ({e^, e^, e z } is the canonical 
basis) 

d - d = -M s e z < } . (5) 
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IV. DYNAMICS OF THE CONDUCTION ELECTRONS 

A detailed derivation of the balance equation for the spin density Sj(x,t) is a many particle problem [24]. We use 
the phcnomcnological expression for the spin current density [18] . 

j fc (x,t) = -/w fc (x,t)E(*)-D(x)V<fo*(x,i) . (6) 

j fe is the spin current density for electrons with spin-polarization along the k— axis, fj, is the electron mobility, 
which we assume as material-independent, and E(t) is the time-dependent electric field. -D(x) stands for the 
material-dependent diffusion constant and #Sfc(x, t) = Sfe(x, t) — s e k q (x,t) is the non-equilibrium spin density (spin- 
accumulation), s e? (x, t) is the space- and time-dependent equilibrium spin density. We compute the latter using 
the SDDE, as explained in the next section. Eq. Q is in general valid for ferromagnetic as well as for nonmagnetic 
materials. Because of V • E = inside the metal, we obtain the spin drift-diffusion equation (SDDE) [T51 |2"5] . 

8s lei 

— = — s x B + (VD • V)<5s + DA5s + 
at m 

^ (E - v)s+ fS) s/ - (7) 



|e| is the elementary charge and m the electron mass (/xo is the permeability of the vacuum). For the spin flip term 
we make a spin-relaxation-time ansatz, 

dtj sf r(x) 

with the space dependent relaxation time r(x). For simplicity we assume an isotropic r inside each layer. In Eq. ufh 
the magnetic induction is related to the magnetization and an external field trough B(x, t) = /io(H(x, t) + M(x, tj). 

A. Spin and charge currents 

From now on we will consider quasi-one-dimensional systems along the x-axis, such as sketched in Fig. [T] Using 
Eq. (|| we obtain for the spin current I s 

I s (x,t) = -A^E(t)s(x,t) + D(x)^(x,t^j . (9) 

Here A is the cross-section of the sample, and E(f) = E(t)e x . In the drift-diffusion model the charge current density 
is given by [2"rJ] 

j(x,t) = n\e\fxE(t)-D(x)— . (10) 

We assume homogeneity, = to obtain 

j(t) = n\e\fiE(t) . (11) 

A related quantity is the drift-velocity Vd, used for numerical computations, defined by t>d(t) = — j(t)/n\e\. For 
electrons j and Vd have opposite sign. Vd > 0, j < means electrons (spin-carrier) move in the positive x— direction. 

B. Equilibrium spin density 

Because we consider an arbitrary movement of the analyzer magnetization vector m(t) we have a time-dependent 
equilibrium spin density (we neglect spin pumping processes induced by moving magnetization [171 127) ). For a fixed 
time t the equilibrium spin density s eq (x,t) is space-dependent, with a first order approximation (isolated layers), 

s eq (x,t) =8= { PWT' ' (12) 

V ' (i) I 0, x £ N. y J 
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This expression reflects our choice of dimensionless spin density s in Eq. ^ and Eq. Q. To obtain the equilibrium 
spin density in the complete structure we use the general stationary solution Eq. ( |14[ ) presented in the next section, in 
which we use the first order expression Eq. ( 12 ) for s eq (x, t). We use the boundary conditions for transparent interfaces 
in F/N-junctions [55]: s(x,t) and l s (x,t) are continuous. For E(t) = we obtain in second order s eq (x,t) |/ 2 ). In the 
following we omit the subscript |( 2 )- Note that, for non-collinear magnetic layers, all components of s eq (x,t) ^ in 
general. 



C. Stationary solution of the SDDE for constant parameters 



For a given layer, we consider Eq. Q for constant current, constant material parameters and time- and space- 
independent magnetic field. We use s eq = s given by Eq. Q. Setting ds/dt = in Eq. 0, we have for a one 
dimensional structure the equation 



D a"(x) + fiE s'(x) + uj s(x) x b x 



s(x) 



. 



(13) 



Here we have defined bi = B/|B|, and B = (m/|e|)wbi, with ui the Larmor frequency. The general solution of 
Eq. (13) is a quite lengthy expression, containing 6 integration constants, denoted as c% . . . c%. To find it we split s(x) 
into two parts, one part parallel to the magnetic field, and the other perpendicular to it, 



s(x) = S II (x) + s±(x) 



(14) 



We define an orthonormal, positive oriented basis {bi, b 2 , b 3 }. One finds for the parallel part (where Id = I^Et is the 
drift length with sign determined by E), 



Sii (x) = bi < ci exp 



x[l d + y/l* + 4A 2 ] 



2A 2 



c 2 exp 



2A 2 



s . 



(15) 



By (a;) does not depend either on |B(t)| or the saturation magnetization. The second part is given by 

_(x) = h 2 \c ?> G A {x) + c 4 G 3 (x) + c 5 G 2 (x) + 
c 6 Gi{x)} + b 3 [c 3 G 3 {x) - c 4 G 4 (a;) - 
c 5 Gi(x) + c 6 G 2 (x)} . 



(16) 



Here the functions Gi(x), i — 1, ...4 are given in Appendix |A| They depend on the magnetic field and the electric 
current, not indicated here to simplify the notation. Eq. ( 14 ) with Eq. (151 and Eq. ( 16 ) present the complete solution 
of Eq. ( 13 1 used in our numerical simulations. We make the following remarks: 

(i) The solution of of the SDDE for spin-orientation-dependent material parameters is straightforward. 

(ii) Using this solution one can study different boundary conditions when linking layers. 

(iii) Because the solutions for spin densities parallel and normal (to the magnetic field) can be separated, it is 
immediately possible to refine the model using different times T\ and r 2 for spin relaxation and dephasing. 



D. Validity of the quasi— static solution 



Here we develop a scheme to estimate the errors from our quasi-static approach. We use the stationary solution 
from the previous section to compute the spin density for a time-dependent current and magnetization vector. In 
general this approximation is valid as long as the variation of j(t) and m(t) is slow compared to the shortest relaxation 
time t (quasi-static time-evolution, QSE). A more rigorous estimate of the accuracy of the QSE in comparison with 
the solution of the full time-dependent equation is a non-trivial task. This is due the different relevant processes and 
time-scales in the different layers. To get a quantitative picture we set s(x,t) = s qs (x,t) + 5s qs (x,t), where s qs (x,t) 



FIG. 1: (color online) Geometry of the seven-layer system. The outside layers act as spin-carrier (electron) reservoirs with 
polarization V = 0. The regions Pi and P2 are the two polarizers, and A is the analyzer layer whose magnetization is to be 
manipulated. 



denotes the quasi-static solution Eq. (14 1 and Ss qs (x,t) the deviation from the exact solution, denoted as s(x,t). For 



Ss qs (x,t) we have inside a single layer the equation 



^ = Mfc,. x B 
at m 



M(E • V)5s qs - 



DAds n 



§S n 



- s„ 



(17) 



The inhomogeneity is denned as 



ds qs dj f dm 



dj dt 



dt 



(18) 



and is the source for a non- vanishing Ss qs . Let us discuss the spin-relaxation in the ferromagnetic layers. Here the 
typical relaxation time is r ps 1 ps and it is reasonable to neglect, in a first approximation, the Larmor, diffusion, and 
drift terms. The Larmor term is of the order of 1/w, the diffusion is characterized by a time scale = I 2 /D = t(1/X) 2 , 
where I is a characteristic finite length (layer thickness). The drift term goes as Tj — l/\vd\- For I = 3 nm (analyzer 
thickness as a worst case) this gives ps 0.3t and Tj ps 0.03 ns for j ps 10 s A/cm 2 . If we integrate Eq. (171 under 
this assumptions we find as a first-order correction (for a ferromagnetic layer), 



Ss q 1 J(x,t) = - / dt'e-^/ T s st (x,t') 
Jo 



(19) 



If we consider now a spacer layer (x± . . .x-% in Fig. [TJ and compare r with = t(1/X) 2 using the parameters given in 
Tab. I we can see that <C t. Diffusion is dominant in the spacer-layers. It occurs on a time-scale ps 10~ 3 ps. In 
fact, this quite different time-scales in different layers are the reason why an integration of Eq. Q by a discretization- 
procedure used in usual PDE-toolboxes leads to numerical problems. To obtain an estimate of 8s q 1 s \x,t) inside the 



spacer-layer we solve Eq. (17 1 with boundary-conditions given by Eq. (19). Numerical results of this strategy to 



estimate the accuracy of the QSE will be given below. 



V. SEVEN-LAYER SYSTEM 



Fig. [JJ shows the seven-layer structure, for which we apply the general formalism. We select this system because 
such structures where used for low-critical current experiments |14| 115) . This allows testing of the present approach. 
As indicated in the figure we use two opposite aligned polarizer-layers, polarizer Pi points in the +z and polarizer P 2 
in the —z direction. P\ defines the parallel position of the analyzer A, where a small deviation from P\ is needed for a 
non-vanishing initial-spin torque. Both polarizers have the same material and geometric properties. As a consequence 
the dipole field Eq. |5]) vanishes exactly at the position of the analyzer A and the anisotropy field Eq. ^ produces 
two energetically equivalent stable positions (degenerate two-level system). We expect and proof that, if a current 
pulse switches the magnetization from P— >AP, then —j(t) does the inverse operation, AP— s-P. 



G 



mat. 


A [nm] 


T [lis] 


M 3 [A/m] 


uj [GHz] 


V 


Cu 


450 


0.024 











Fe 


5 


0.001 


17 x 10 5 


230 


0.45 


Py 


5 


0.001 


8 x 10 5 


110 


0.37 



TABLE I: Material-parameters used in the simulation. 



A. Numerical strategy 



All computations are done with the help of Mathematica. We use the solution Eq.(14|-(16) to compute the time 



and space-dependent spin density inside of each layer for given direction of B and current j . The solution of the total 
system requires the determination of all integration constants C\ . . . C36 . Whereas the boundary conditions (continuous 
spin density and spin current density) are formulated analytically, the solution is computed numerically as a function 
of m and j. The spin current density and the spin torque in Eq. p| are then calculated self- consistently using 
The last step requires the numerical solution of the LLGE, Eq. p|. 



B. Optimized switching procedure 

We now address the switching of the analyzer magnetization (for optimized switching using external magnetic fields 
see [H]). We first note that, due to the non-linearity in m of the LLGE, it is impossible to identify a single current 
pulse profile which switches both from P— > AP and AP— >■ P (initial-state-independent switching). However, using 
the symmetry of the structure, one can identify a current pulse profile which, when changing the current direction 
only, promotes both processes. 

To find a simple pulse shape which performs the desired task it is convenient to use an optimization procedure based 
on a suitably defined cost-functional J [3D]. We set 

J = ||m(tj) — my 1 1 , 0<J<2, (20) 

where my is the target magnetization and m(i/) is its actual value at the prescribed target time tf. We choose the 
time-dependent current as 

j(t;Xi,X 2 ,X 3 ) = 

X A expf.-^{t-t f /2) 2 ^+^X lB m(lirt/t f ) , (21) 

with 3 variational parameters X±, . . . X3 (one can use also more parameters. Global optimization algorithms however 
work best with a few parameters). The Gauss-pulse, characterized by Xa, Xb, is selected by hand such that it is 
sufficient to switch the magnetization from P— >AP. It is used as a reference pulse. However, to steer the magnetization 



in the prescribed time additional current contributions are needed. The additional terms in Eq. (21) are constructed 
to ensure that at the end points of the control-time interval [0, t/] the current vanishes for arbitrary Xi.,2,3. To find 
the minimum of J(X\, . . . X3) a standard line search method or genetic algorithm can be used. 



C. System— Parameters 



We use the material-parameters typical for a Cuqo /Fei5/Cu3/Py2/Cu3/Fei5/Cu 00 (in nm) multilayer-system. The 
relevant material-parameters are listed in Tab. I [HI [31] [32] • We use a material-independent electrical conductivity 
and free-electron-density of n — 84/nm 3 . We obtain a microscopic expression for the coupling constant £ given by 
£ = ~~ 2mM d ~ ~0.969/d, where d — X3 — X2 is the analyzer thickness. For the Gilbert damping parameter we set 
a = 0.01 33 . The direction of the anisotropy field is chosen as n = (0, sin(y), — cos(<p)), with ip = 0.97T, and its 
modulus co an — 2 GHz. 
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FIG. 2: (color online) a) Electric current and b) 2-component of the magnetization in the analyzer versus time. Associated 
quantities are plotted in the same line style. For decreasing current the switching time increases. The electric current is plotted 
as — j according to a positive drift velocity. 



D. Results 



1. Switching into constant current 

For a first example we consider the dynamics of the seven-layer system in Fig. [T] for a current that we switch on 
according to j(t) — j'q(1 — e~'/ T ), with T = 0.5 ns. We have integrated the LLGE for different values jo, as shown in 
part a) of Fig. [2] Part b) shows the z— component of the magnetization as a function of time. In all three cases the 
magnetization switches from P— >AP, however, the lowest current leads to a switching time of more than 100 ns. These 
investigations agree well with basic experimental results in the literature [21 [TS] : the critical current \j c | is of the order 
\j c \ ~ 10 6 A/cm 2 , for the parameters chosen here, and depends on the saturation magnetization, Gilbert damping, and 
anisotropy field |16j . As seen in Fig. [2j switching into a constant spin-polarized current leads to damped oscillations 
of the magnetization vector. Above the critical current, they result in a flipping of the magnetization vector into the 
new (AP) equilibrium position. For currents \jo\ < \j c \ one induces damped oscillations without switching. We should 
remark that the equilibrium-positions of m for a constant (spin-) current are no more given by the directions of ±n, 
but there is small deviation due to the spin current, however, not resolved in Fig. [2} 

The seven-layer structure with antiparallel polarizer-orientations is crucial for the occurrence of low \j c \. Com- 
putations for parallel polarizer orientations (-P1II-P2) give vastly different critical currents for P— !• AP and AP— >P 
flips. Fig. [3] reveals the reason for this result. For anti-parallel orientation of the polarizers the z— component of the 
spin density shows a large gradient inside the analyzer layer. As a consequence large spin currents can be generated 
compared to parallel oriented polarizers. In fact for a simplified model with vanishing dipole field (for sample radius 
R — » 00) and parallel polarizers the critical current is \j c \ > 10 s A/cm 2 for this structure. As investigated, switching 
times for the analyzer magnetization tend to decrease with increasing |jo| |34j . 



2. Optimal pulse-sequences 



We now consider the problem of switching of the magnetization m using an optimized time-dependent electric 
current, where we set the switching time to tf = 5 ns. The first current pulse should switch the magnetization from 




x (nm) 



FIG. 3: (color online) Equilibrium spin density for two different polarizer alignments. The blue-dashed line corresponds to 
parallel orientations of the two polarizers, the black-solid line is for anti-parallel orientation, as used in low-current spin-torque 
experiments. The analyzer is at the position m = n and therefore small perpendicular components of the spin density are 
present. 



P— ^AP. Initial and desired final value of the analyzer magnetization m(i), respectively, are 

m(0) = n and m(t/) = mj = — n . (22) 



A numerical minimization of Eq. (20), limiting ourselves to the pulse shape Eq. 21 gives as a result the first pulse 
shown in Fig. [5] We stopped the computation when the cost functional was J s» 0.006. This means that the optimal 
control pulse, rather than relying on intrinsic Gilbert damping, actively drives the magnetization precisely into the 
target state AP; likewise for the back flip, see Fig. [5](c). Note that the pulse shape is chosen such that the current is 
zero at the boundaries of the time interval. To ensure that the magnetization remains in the AP state after the first 
flip, a few ns later we apply the same pulse once more. Only a weak deviation from the equilibrium position in form 
of a few damped oscillations are visible demonstrating stability, see Fig. [5] However when we apply the same pulse 
profile with opposite current direction we switch the magnetization back from AP— >P. In addition to m(t) we have 
plotted in Fig. [5] the time-dependent spin density during the first current pulse. The rows (a) and (b), respectively, 
show the equilibrium spin density and its deviation from equilibrium inside the multilayer device. One observes the 
degree to which the equilibrium spin density depends on the time-dependent magnetization m(t): due to the choice 
of the magnetization of Pi and P2 (as collinear) only the z-component of the spin density shows significant deviation 
from equilibrium. The order of magnitude of the deviation of x- and y-components is of the order of the error made 
by the QSE. The non-equilibrium spin density as function of time is influenced by the actual position of m(£) and 
the current j(t), as well as the magnetization of Pi and P2 ■ 



3. Error estimate 



We have used the results from the previous section to test the numerical validity of the stationary solution as 
discussed in Sec. IV D Fig. [6] shows the estimate for the deviation of the z— component of the spin density in selected 
parts of the structure. The solid line is for the center of Pi , while the dashed line is for the center of the spacer layer to 
the left of the analyzer. The figure shows that [£s^] z , depending on position, is of the order of 10~ 5 — 10~ 4 , compared 
with s e z q « 0.1-0.3 (see Fig. [3). The dominant contribution in Eq. Q comes from the moving magnetization, whereas 
the current contribution is neglibile. For the other components we obtained similar results regarding relative errors. 



VI. CONCLUSIONS AND OUTLOOK 



We have presented a self-consistent model for magnetization switching by spin-polarized electric current in metallic 
ferromagnetic heterostructures. Our method is founded upon an analytic solution of the stationary spin drift-diffusion 
equation (SDDE) for each layer using constant material parameters, electric current, and magnetic field. Matching 
layers, using continuity of spin density and spin current density at the interfaces as boundary conditions, we obtain 
an analytic solution for the spin density of the entire heterostructure. Making a quasi-static approximation in which 
the time dependence of the spin density depends on time solely via the electric current and net magnetic field, the 
time evolution of the spin density is computed in parallel to the Landau-Lifshitz-Gilbert (LLGE) equation. Both 
equations couple via the spin torque effect and the time-dependent magnetization in the SDDE. This method al- 
lows for an efficient and robust mathematical description of the coupled carrier spin and magnetization dynamics in 
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a) c) 

-j x 10 7 (A/cm 2 ) 




FIG. 4: (color online) a) Optimized time-dependent electric current for switching the analyzer magnetization from P— >AP 
and vice versa. Three 5 ns pulses with the same shape are applied. The first switches from P— >AP, the second is used to test 
stability, and the third switches back AP— >P. b) ^-component of the analyzer magnetization vector m(t) as a function of time, 
c) Plot of the three-dimensional trajectory of m(t) during the first (black, solid) and last pulse (blue, dashed). 



metal/ferromagnet heterostructures. Because the model is based on a completely analytic solution of the stationary 
SDDE for given electric current and magnetic field for each layer, it is applicable to heterostructures of high complex- 
ity, for example for tilted polarizers or structures exposed to external magnetic fields [35^ . 

We have demonstrated the efficiency of this semi-analytic approach by investigating a seven-layer system with an- 
tiparallel oriented polarizers, as studied in recent experiments, and computed optimized current pulses to switch the 
magnetization from P— >AP— >P in specified time of 5 ns. As expected for the system under investigation, the obtained 
current densities are in the range of 10 s A/cm 2 , with a critical current of about 10 6 A/cm 2 . Using optimal control 
theory, we identify solutions for current profiles which allow for precise switching in predetermined switching times. 
We provide and discuss one example. 

Furthermore, a detailed investigation of the validity of the quasi-static time evolution of the SDDE is given. It 
confirms excellent accuracy for the example of the simulated seven-layer heterostructure. 

Several future applications of the presented formalism can be envisioned. A combined variation of material- and ge- 
ometric parameters to obtain optimal current pulses with low critical currents. A description of thermal fluctuations 
using temperature-dependent effective (Langevin-) fields in the LLGE (via the spin torque in the SDDE) and the 
search of "thermally robust" current pulses by averaging over many field configurations. 
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FIG. 5: (color online) a) The first row shows the equilibrium spin density for the switching process P— > AP for the first pulse 
in Fig. [4] It depends on m(t). b) Non-equilibrium spin density induced by the electrical-current pulse. The computed values 
for 5s x , Ss y , however, are at the limit of the accuracy of the QSE. 
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FIG. 6: (color online) Numerical estimate of the error within the QSE relative to an exact treatment of the SDDE. The inset 
shows the locations where we compute Ss 
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Inside the polarizer (red-solid line) we use Eq. ( 19 1 to estimate the deviation 
from the exact result, whereas inside the spacer-layer, Eq. (|17[) is integrated numerically. 
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Appendix A: Stationary solution of the SDDE 



Here we summarize the remaining analytic expressions for the stationary solution and constant material parameters 
as presented in Sec. IV C We use the dimcnsionlcss quantities k := lot and p := ld/X. Further we define 



- 4in} 



4 + p 2 



\ 8 N 



4 + p ; 
4k 



(Al) 



Using the auxiliary functions, 



b = Im[\/4 + p 2 +Mk] 



\ 



4 + p 2 



\ 



4 + p 2X ' 



4k 



Fj (./ ) = f ~* cos ( ^ ) siiib 



V2aJ 



(A2) 



(A3) 



F 2 (x) = e ^ a cosh (^yj sin 



(A4) 



F 3 (a;) = e 2A cos 







) C0S11 (2X) 


V2Ay 





(A5) 



bx \ / \ 



the four dimensionless functions Gi(x), entering in Eq. (16) are: 



(A6) 



Gi(x) = [-4a 3 pn + 4ap K (A + 3b 2 + p 2 )] Fi(x) + 
[4b 3 pn + Abpn{A - 3a 2 + p 2 )} F 2 (x) - 
4 K (a 2 + 6 2 )(-4 + a 2 - b 2 - p 2 )F 3 (x) - 
8ab(a 2 +b 2 )nF 4 (x) , 



(A7) 



G 2 {x) = ap [a 4 + (4 + 5b 2 + j o 2 )(4 - 2a 2 + b 2 + p 2 )] F t {x) + 
bp [5a 4 + (4 + b 2 + p 2 ) 2 - 2a 2 (12 + 56 2 + 3p 2 )] F 2 (x) - 
(a 2 - b 2 ) [(a 2 + b 2 -p 2 - 4) 2 - 4aV] F 3 (x) + 
4a6(a 2 - b 2 )n [4 - a 2 + b 2 + p 2 ] F 4 (x) , 



(A8) 



G 3 (x) = [2ia 2 bn - 8bn{4 + b 2 + p 2 )} F x {x) + 
[24ab 2 K + 8ok(4 - a 2 + p 2 )] F 2 (x) , 



(A9) 



Gi{x) = [-8a 3 k + 8a K (4 + 3o 2 + p 2 )} F t (x) + 
[% 3 k + 86k(4 - 3a 2 + p 2 )} F 2 {x) . 



(AlO) 
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The integration of the normal component of Eq. ( 13 ) requires the solution of two second order differential equations. 
It is advantageous to transform this two second order equations into four first order equations and solve this system 
by matrix exponentiation. This procedure, after some simplifications, leads to the four functions Gi(x), which build 
the fundamental solution. 
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